Comparison of Density Functional Approximations 
and the Finite-temperature Hartree-Fock Approximation in Warm Dense Lithium 



Valentin V. KarasievQ Travis Sjostrom, and S.B. Trickey 
Quantum Theory Project, Departments of Physics and of Chemistry, 
P.O. Box 118435, University of Florida, Gainesville FL 32611-8435 
(Dated: 19 July 2012) 

We compare the behavior of the finite-temperature Hartree-Fock model with that of thermal den- 
sity functional theory using both ground-state and temperature-dependent approximate exchange 
functionals. The test system is bcc Li in the temperature-density regime of warm dense matter 
(WDM). In this exchange-only case, there are significant qualitative differences in results from the 
three approaches. Those differences may be important for Born-Oppenheimer molecular dynamics 
studies of WDM with ground-state approximate density functionals and thermal occupancies. Such 
calculations require reliable regularized potentials over a demanding range of temperatures and den- 
sities. By comparison of pseudopotential and all-electron results at T = OK for small Li clusters of 
local bcc symmetry and bond-lengths equivalent to high density bulk Li, we determine the density 
ranges for which standard projector augmented wave (PAW) and norm-conserving pseudopotentials 
are reliable. Then we construct and use all-electron PAW data sets with a small cutoff radius which 
are valid for lithium densities up to at least 80 g/cm 3 . 



I. INTRODUCTION 



Warm dense matter (WDM) encompasses the region 
between conventional condensed matter and plasmas. 
WDM occurs on the pathway to inertial confinement fu- 
sion and is thought to play a significant role in the struc- 
ture of the interior of giant planets. The theoretical and 
computational description of WDM is important for un- 
derstanding and performing experiments in which WDM 
is created [l( . Two parameter ranges which are very dif- 
ferent from those in standard condensed matter physics 
characterize WDM: elevated temperature (from one to a 
few tens of eV) and high pressure (up to thousands of 
GPa). These ranges are challenging computationally be- 
cause the standard solid state physics methods become 
very expensive (due to high temperature) or standard ap- 
proximations used in those methods cease to work (due 
to high material density). From the plasma side, the 
temperature and pressure are not high enough to employ 
classical approaches. 

A combination of a quantum statistical mechanical de- 
scription of the electrons, and classical molecular dynam- 
ics for ions is a standard theoretical and computational 
approach to WDM at present. Usually the quantum sta- 
tistical mechanics is handled via finite-temperature den- 
sity functional theory (ftDFT) There is a sub- 
stantial literature, too large to review here, about such 
calculations at zero temperature via Born-Oppenheimer 
molecular dynamics (BOMD) or Car-Parrincllo MD, 
with DFT implemented via the Kohn-Sham (KS) pro- 
cedure for the electronic degrees of freedom. The per- 
tinent point is that the same techniques can be applied 



to the finite-temperature case The combination, 

called ab initio molecular dynamics (MD), is computa- 
tionally costly at high temperature (for a given density) 
because of the large number of partially occupied KS or- 
bitals which must be taken into account. 

The great majority of the reported finite-temperature 
ab initio MD calculations use zero-temperature 
exchange-correlation (XC) functionals, E xc , with 
Fermi thermal occupations to construct the electron 
density. In such calculations, the only T-dcpendcncc in 
the XC contribution to the free energy J-" xc is through 
the T-dependence of the electron density: 



Jxc[n(r,T),Tl«B»cKr,T)] 



(1) 



with n(r, T) the electron number density at temperature 
T. 

Most ftDFT calculations with ground state XC func- 
tionals seem to have been done with the Vasp [2(| or 
Abinit [2l| codes using either the local density approx- 
imation (LDA) for £ xc |2a - l24| or the Perdew-Burkc- 
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Ernzerhof generalized gradient approximation (GGA) 
functional [2a ]. 

The orbital-free density functional theory (OF-DFT) 
treatment of electronic degrees of freedom is a less ex- 
pensive alternative to orbital-dependent methods such as 
KS. OF-DFT in principle provides the same quantum- 
mechanical treatment of electrons as KS DFT, but the 
lack of accurate orbital-free approximations for the ki- 
netic energy functionals has limited the use OF-DFT, 
even at standard conditions. In contrast, the high density 
of the WDM regime is favorable for use of the OF-DFT 
approach, which is a motive for developing functionals. 
The standard KS approach clearly must be used to test 
and calibrate such OF-DFT functionals. The limitations 
and consequences of various choices in those thermal KS 
calculations have not seen much detailed attention how- 
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ever. Two closely related sets of potentially significant 
issues occur. 

First, the use of ground-state functionals in a ftDFT 
calculation inevitably raises a topic for fundamental 
DFT, namely, the adequacy, accuracy, and scope of Eq. 
([1]). Relative to the number of calculations, there are 
comparatively few studies to assess this approach against 



comparatively tew studies to assess tnis approacn against 
others @, !, 0, BM HI • Ref. shows that the maxi- 
mum density of the Al shock Hugoniot is increased about 
5% or less by use of a temperature-dependent functional 
of the Singwi-Tosi-Land-Sjolander (STLS) type |27| . Ref. 
l2ril made essentially the same comparison but with re- 
spect to simple Slater exchange (in Hartree atomic units) 
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and with the added complication [for the purpose of as- 
sessing Eq. ([1])] of use of an OF-DFT approximation. 
Ref. [8| compared calculations for ground-state LDA and 
PW91 GGA [11] functionals. Faussurier et al. [3] com- 
pared the electrical conductivity of Al computed with 
the T-dependence from classical-map hypernetted chain 
scheme [29[ versus ground-state LDA. They concluded 
that the effects on conductivity are small in the WDM 
regime but become increasingly important as the energy 
density increases. Wiinsch et aZ.[lj| reversed the perspec- 
tive and used ftDFT calculations with a ground-state 
XC functional to calibrate hypernetted chain approxi- 
mations, hence assumed the validity of Eq. (fTJ). Vinko 
et al. compared ground-state GGA calculations of free- 
free opacity for Al with an RPA model and found semi- 
quantitative agreement at lower photon energies with in- 
creasing disagreement at higher ones, all over the range 
< T < 10 eV [l7| . As an aside, we note that the same 
issues of use of ground-state approximate XC function- 
als in a T- dependent context can arise in average-atom 
models ^Q^. 



The second set of issues involves computational tech- 
nique. The primary focus is control of the effects of pseu- 
dopotentials (or regularization of the nuclear-electron in- 
teraction). These are ubiquitous in the highly rehned 
codes in use for both WDM and ground state calcula- 
tions. Clear insight into the behavior and limitations of 
functionals requires that the regularized potentials not 
introduce artifacts of their own. The challenge is to test 
those potentials against high-quality all-electron (AE) re- 
sults over the appropriate density range. 

An obvious issue associated with pseuodpotentials is 
the effect of a finite core radius upon compressibility 
(hence, equation of state). Ref. |35| shows that a norm- 
conserving pseudopotcntial for boron with the standard 



cutoff radius (r c = 1.7 Bohr) is not transferable to the 
high material density regime. In that work, the au- 
thors built an "all-electron" pseudopotential with small 
r c = 0.5 Bohr and tested its transferability to very high 
material density by comparison with the Thomas-Fermi 
(TF) limit calculated using an average-atom model ■ 
Another issue is the extent to which removal of core 
electrons has an unphysical effect on the distribution 
of ionization. A related issue is the effect that remov- 
ing core levels has on Fermi-Dirac occupation numbers. 
At fixed density, such core levels should be progressively 
depopulated with increasing temperature. Does the de- 
population of pseudo-density levels behave correctly? A 
significant computational practice issue is the minimum 
magnitude threshold for retention of occupation num- 
bers. That threshold is directly related to basis set size 
or, equivalently, the plane- wave cutoff. We know of only 
one study of any of these questions [36| . In it, all-electron 
calculations with the full-potential linearized muffin-tin 
orbital methodology were used to benchmark projector 
augmented wave (PAW) calculations with a plane wave 
basis. Two metals, Al and W, were treated at T ^ OK. At 
least for W, it appears that different XC functionals were 
used for the comparison. Additionally, Ref. [36| used the 
free-electron expression for the non-interacting electronic 
entropy, rather than the proper explicit dependence on 
occupation numbers ji\ 

S s = -k B Y^{fi In fi + (1 - /*) ln(l -/«)}. (3) 

i 

Despite these differences, to the extent that their topics 
and ours overlap, the findings are consistent. 

To establish a basis for comparison, first we consider 
the issues of regularized potentials. We consider both or- 
dinary pseudopotentials (PPs) and the pseudopotential- 
like PAW technique. Those tests are against all-electron 
(bare Coulomb nuclei potential) calculations for small Li 
clusters of bec symmetry. We establish a PAW which 
demonstrably is reliable for the density range of inter- 
est. Then we study the behavior and limits of the 
use of ground-state X functionals in ftDFT by compari- 
son of finite-temperature Hartree-Fock (ftHF) and DFT 
exchange- only results. For clarity of interpretation, all 
the bulk solid calculations reported here were performed 
at fixed ionic positions corresponding to an ideal bec 
structure for Li. 



II. CODES 



We used the atompaw code 371 to form the PAWs 



For periodic systems, we used three codes, Abinit vers. 
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Vasp vers. 5.2 [2Cj, and Quantum-Espresso 
381 ] . All three are plane- wave, PP codes. All 
three also implement PAWs. Abinit and Quantum- 
Espresso are open source. Technical details of the ftHF 
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calculations arc discussed below. For the all-electron cal- 
culations on finite clusters, we used conventional molec- 
ular gaussian basis techniques as embodied in the Gaus- 
sian 03 program J^. 



III. REGULARIZED POTENTIALS 

Diverse PP techniques commonly are used in KS cal- 
culations to reduce computational cost by excluding the 
core electrons from the self-consistent field (SCF) pro- 
cedure and to regularize the singular external potential 
in order to use an efficient, compact plane wave basis 
set. Excluding core electrons implicitly invokes the frozen 
core approximation (i.e., the omission of core electrons 
from the SCF procedure) . That approximation generally 
is well-justified in standard conditions. There, the core 
electrons are uninvolved in chemical bonding and their 
state is essentially independent of the chemical environ- 
ment. The validity of this justification is not obvious for 
the WDM regime. In it, all electrons become important 
for correct evaluation of the Fermi occupancy at high 
temperature and correct description of the electron den- 
sity at high external pressure. As a consequence, it is 
mandatory to include at least some core electrons in the 
solution of the relevant Euler equation (DFT or finite- 
temperature HF) in the WDM regime. For light atoms 
this may mean an all-electron PP. Those are, of course, 
a particular form of regularized potential. 

Generation of PPs usually is characterized by cutoff (or 
pscudization) radii, r c . Values of r c are a compromise be- 
tween softness of the PP (for compactness of plane wave 
basis sets) and correct description of the one-electron 
orbitals close to the nucleus. Standard PPs are devel- 
oped for use under near-equilibrium condensed matter 
and molecular conditions, hence their transferability to 
the WDM regime needs to be explored. For example, 
commonly r c is assumed to be somewhat smaller than 
half the nearest-neighbor distance between atoms so that 
there is no core overlap. There is no guarantee that such 
equilibrium prescriptions are satisfactory for WDM stud- 
ies. 



A. Basic PAW formalism 

PAW concepts are summarized in Ref. 0. We out- 
line the relevant points here. The PAW valence electron 
energy is comprised of a pseudo-energy evaluated using 
a smooth pseudo-density and pseudo-orbitals plus atom- 
centered corrections. An energy correction centered on 
atom a is evaluated using an augmentation sphere of 
radius r°. Within each sphere, the correction replaces 
the valence pseudo-energy of atom a, E%, by the valence 
energy E% generated from the valence part of the all- 



electron atomic density 



En — En 



E 



(4) 



Detailed descriptions_of each term in Eq. (jU are given, 

Here the issue is treatment of 
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for example, in Ref. 
core density contributions to the XC energy, as discussed 
in that reference. In the scheme due to Blochl [42| , the 
XC energy is expressed as 



EL 



(5) 

where n a and n° are atom-centered valence and core elec- 
tron charge densities corresponding to all-electron atomic 
orbitals, n a and n° are atom-centered valence and core 
electron pseudo-densities, and h, h c are total valence and 
core electron pseudo-densities. The idea behind Eq. ([5]) is 
that the third term, which corresponds to atom-centered 
contributions of pseudo-densities (evaluated within aug- 
mentation spheres, radii r"), cancels the corresponding 
atom-centered pseudo-density contributions (evaluated 
over all space) in the first term, and the canceled contri- 
bution is replaced by the second term, which is evaluated 
with atom-centered all-electron densities (again within 
the augmentation spheres only). 

The Kresse scheme [43[ introduces a valence compen- 
sation charge density, n, as well. Its purpose is to repro- 
duce the multipole moments of the all-electron charge 
density outside the augmentation spheres [41|. For the 
XC contribution, n is added to the pseudo-densities in 
the functionals in Eq. (O to give 



E xc — E 



-n a c }-E xc [n a +n a c - 



(6) 

This procedure can cause problems with GGA XC func- 
tionals; see Ref. H^. 

There are what are called all-electron PAWs, which 
in essence are regularized potentials for all-electron cal- 
culations. In customary notation, an "N-electron" 
PAW retains N electrons in the valence. Thus a 3- 
electron ( "3e~" ) PAW calculation for Li is an all-electron, 
regularized-potential calculation. 



B. PAW and high density lithium 

We tested the PAW approach by calculating the pres- 
sure of bec Li over a large range of material densi- 
ties, from approximate equilibrium, pu = 0.5 g/cm 3 , to 
PLi = 25.0 g/cm 3 , all at T = 100K. (The equilibrium 
density from simple Slater LDA all-electron calculations 
is 0.54 g/cm 3 , or lattice constant 6.59 Bohr, close to the 



experimental value; see Ref. |44j. Newer LDAs give some- 
what contracted results; sec below.) Three different PAW 
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data sets were used for each LDA and GGA exchange- 
correlation functional: (i) the standard set with compen- 
sation charge density included from Ref. Hfl, (ii) a set 
with the same cutoff radius (r c = 1.61 Bohr) but with- 
out compensation charge density, and (iii) a set we gen- 
erated with r c = 0.80 Bohr and no compensation charge 
density. The Perdew-Wang (PW) and Perdew-Zunger 
(PZ) LDAs [H, HI and Perdew-Burke-Ernzcrhof GGA 
(PBE) XC functionals were used. 

The upper segment of Table [J compares the calculated 
bec Li equilibrium lattice constants and bulk moduli for 
the various combinations. These were done with Abinit. 
The lattice constant and bulk modulus were obtained 
by fitting the calculated total energies per cell to the 
stabilized jellium model equation of state (SJEOS) form 
[47I ]. One sees that the exclusion of the compensation 
density slightly decreases the lattice constant for both 
PW and PBE functionals. The results are essentially 
unchanged when the r c value is decreased to 0.80 Bohr. 

Table fl] also summarizes results obtained using both 
Quantum-Espresso and Vasp. The lattice constant 
and bulk modulus again were obtained via fitting to the 
SJEOS form in all cases. The results for Vasp come from 
using the PAW pscudopotcntials supplied with the code 
itself. There is excellent agreement between Quantum- 
Espresso and Abinit results when the same 3e~ PAW 
data set is used. The Vasp PBE 3e~ results do not 
agree as well, consistent with the findings of Ref. re- 
garding the effects of the valence compensation charge 
density contribution. Two LDA PPs also were used with 
Quantum-Espresso, namely the le~ Von Barth-Car 
and 3e~ norm-conserving pscudopotcntials (both taken 
from the Quantum-Espresso web page). The lattice 
constant corresponding to the first of these PPs is un- 
derestimated as compared to other PZ LDA calcula- 
tions, independent confirmation of the importance of the 
3e~ treatment. For the le~ (Vanderbilt ultrasoft) and 
3e~ (norm-conserving) PBE PPs (again taken from the 
Quantum-Espresso web page), the lattice constant is 
slightly overestimated and the bulk modulus is underes- 
timated by the le~ pscudopotential. The 3e~ results are 
in nearly perfect agreement with the PAW data. 

To assess the PAW method for high material density, 
we compared PAW and true all-electron (bare Coulomb 
potential) results for two small lithium clusters with lo- 
cal bec symmetry; see Fig. [1] The interatomic dis- 
tances in both clusters were set equal to the nearest- 
neighbor distance in bulk bcc-Li for densities in the range 
0-5 < pu < 150 g/cm 3 . The all-electron calculations 
were done with the Gaussian 03 code and two basis sets, 
6-311++G(3df,3pd) and cc-pVTZ. For the LDA calcula- 
tions, we used the Vosko, Wilk, and Nusair parameter- 
ization (VWN) [22J]; it is very close to the PZ parame- 
terization and based on the same data. For the GGA 
functional we used PBE. For high densities, pu > 50 



g/cm 3 , we did additional calculations with cc-pV5Z (8- 
atom cluster) and cc-pVQZ (16-atom cluster) basis sets. 
The PAW calculations were done with the Abinit code. 
In it, the clusters were centered in a large cubic super-cell 
of size L. For the standard PAW, we used L = 15 A with 
an energy cutoff 1000 eV, while L = 12 A with energy 
cutoff 3000 eV was used for the small r c PAW. The PZ 
LDA and PBE GGA functionals were used in these cal- 
culations. Note that the difference in behavior between 
PZ and PW is essentially negligible for the purposes of 
this study. 

Fig. [5] shows all-electron and PAW LDA total ener- 
gies for the two clusters as a function of distance corre- 
sponding to the stated bulk density. Fig. [3] shows the 
corresponding GGA results. The behavior of the two 
clusters is quite similar. For the standard PAW data set 
(labeled "(i)" previously), the total energy starts to de- 
viate from the all-electron (AE) values at approximately 
p dust-criti = g q g / cm 3 For thc s t an dard PAW set with- 
out compensation density [set (ii)], the critical density 
pciust-crit2 j g approximately 25 g/cm 3 . In contrast, thc 
PAW with small r c and no compensation density [set (iii)] 
gives essentially perfect agreement with the AE results 
for the whole density range. For densities up to 30 g/cm 3 , 
two basis sets, 6-311++G(3df,3pd) and cc-pVTZ, give 
essentially the same quality results. At high density (50 
g/cm 3 and up), the cc-pVTZ basis set energies lie above 
the values corresponding to the 6-311++G(3df,3pd) ba- 
sis. For those high densities, AE calculations done with 
the larger cc-pV5Z (8-atom cluster) and cc-pVQZ (16- 
atom cluster) basis sets lower the total energy to the 
6-311+-|-G(3df,3pd) level (16-atom cluster) or slightly 
lower (8-atom cluster). Once again there is essential per- 
fect agreement with the set (iii) PAW plane wave results. 

The corresponding PBE GGA comparison of PAW 
and AE results, Fig. [3l shows that the critical densities 
for each PAW data set are almost identical for the 8- 
atom and 16-atom clusters. For the PAW data set (i), 
pfy-st-criti _ g q g/ / cm 3 is lightly i owcr than for thc 

LDA case. For PAW data set (ii), the critical density 
is essentially the same as for LDA 25 g/cm 3 ). Once 
again, the small r c PAW data set (iii) gives good agree- 
ment with the AE results up to the maximum density 
considered (150 g/cm 3 ). We conclude that PAW data 
set (iii), namely r c = 0.80 Bohr and no compensation 
charge, may serve for making reference KS calculations 
in thc high density regime. 

The remaining validation issue is the effect of PAW 
or PP on the calculated pressure. A study [ll| of the 
EOS for warm, dense LiH found that the 3e~ PAW for 
Li in VASP calculations was necessary for T = 2, 4, and 
6 eV and densities twice that of ambient and greater. 
Fig. |4] shows the bulk bec Li pressure as a function of 
material density at T=100 K calculated using Abinit 
with thc same three PAW data sets as before for both 
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TABLE I: Equilibrium lattice constant for bcc-Li, a (Bohr) and bulk modulus, B (GPa). 
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"Compensation charge density is included. 

'LDA: PZ exchange-correlation, nonlinear core-correction Von 
Barth-Car; GGA: PBE exchange-correlation, nonlinear core- 
correction, Vanderbilt ultrasoft pseudopotentials. 

C PZ and PBE semicore state s in valence Troullier-Martins pseu- 
dopotentials. 
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FIG. 1: Bcc Lis (left panel) and Liie (right panel) clusters 
used to test PAW calculations. 
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FIG. 2: All-electron (VWN correlation) and PAW (PZ cor- 
relation) LDA total energies for the Lis and Lii6 clusters. 
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FIG. 3: All-electron and PAW GGA total energies for the 
Lig and Lii6 clusters. 



the PW LDA and PBE XC functional. (Use of the PZ 
LDA functional gives results indistinguishable from those 
from PW LDA on the scale of the figure.) One sees 
that the standard PAW data set (i) starts to overesti- 
mate the pressure at p^ lk - crltl = 6.0 g/cm 3 for LDA 
and at a slightly lower value for PBE. PAW data set (ii) 
produces results which agree with the reference calcula- 
tions (i.e., those from PAW set (iii)) for densities up to 
pbuik-cnt2 _ 25 q g/ cm 3 Comparison of critical density 
values in the clusters and in bulk shows that p^j lk - cntl 
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FIG. 4: Pressure vs. material density from PAW LDA (PW 
correlation) (left panel) and from PAW GGA (PBE XC) (right 
panel) calculations for bulk bcc-Li. 
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panel) for material densities in the range 0.6 — 10.0 
g/cm 3 obtained from Vasp and Quantum-Espresso 
using standard PPs (the PAW provided with the Vasp 
package and the norm-conserving PP taken from the 
Quantum-Espresso web page), and reference results 
obtained with our PAW data set (iii) (for both the PZ 
LDA and PBE GGA XC functionals.) For LDA, we also 
show the earlier all-electron results by Boettger and Al- 
bers [HI . Observe first that our designation of the PAW 
(iii) as a reference is substantiated by the agreement with 
the all-electron LDA calculation. The Vasp le~ PAW 
LDA results start to deviate from the reference values 
at about 3.0 g/cm 3 . However the 3e _ PAW LDA pres- 
sure from Vasp is in good agreement for densities up 
to 6.0 g/cm 3 . Quantum-Espresso results calculated 
with the le~ PZ LDA pseudopotcntial start to deviate 
from the reference results for density between 2 and 3 
g/cm 3 , whereas the 3e~ potential produces results which 
agree virtually perfectly for the full density range. For 
the GGA case, the right-hand panel of Fig. [5] shows that 
the situation is very similar except that both the le~ and 
3e~ PAW Vasp calculations start to deviate from the ref- 
erence results at almost the same density (ps 3 g/cm 3 ). 
Quantum-Espresso le~ calculations overestimate the 
pressure for pu > 0.8 g/cm 3 . However, the Quantum- 
ESPRESSO 3e~ results are in virtually perfect agreement 
with the reference PAW results for the whole range of 
densities. 



p L , <g'< 



p Li (g/cm") 



FIG. 5: Validation of Vasp le~, 3e~ PAW, Quantum- 
Espresso le~, 3e~ PZ LDA (left panel) and Vasp le~, 3e~ 
PAW, Quantum-Espresso le~, 3e" PBE GGA (right panel) 
pseudopotential calculations: pressure as a function of density 
for bcc-Li calculated at T=100 K. 

is slightly lower than pdust-enti ^ A crude linear extrapo- 
lation of the results from PAW data sets (i) and (ii) gives 
an estimated lower bound for the critical bulk density for 
the reference PAW data set p^-critz tQ bfi 8Q g / cm 3 
Additional tests would be needed to get the actual value 
of p^ lk - cnt3 _ Such a determination is not required for 
the present purposes. 

We observe that for fee Al at T=0 K, Levashov et 
aZ.jlBl found that the standard VASP PAW pressures 
began to deviate materially from all-electron values at 
about a compression of seven. Since it was standard 
VASP, presumably that PAW included charge compen- 
sation, hence their result should correspond to our set (i) 
Li results, those labeled "PAW, r c = 1.61 bohr, c.ch." in 
Fig. [4] It is clear that the deviation they found in fee 
Al is at similar but modestly lower compression than we 
find for bee Li. 

Figure [5] compares pressure versus material density 
(T=100 K) for LDA (left panel) and PBE GGA (right 



IV. FINITE TEMPERATURES 

A. Pseudopotentials and level populations 

In finite-temperature calculations (either KS or HF), 
there is non-zero occupation of one-electron levels which 
correspond to empty levels at T = K (virtual states or 
simply "virtuals"). Satisfaction of some computational 
threshold for the smallest non-negligible occupation num- 
ber requires an increasingly large set of those virtuals to 
be considered with increasing T. Concurrently there is 
depopulation of levels fully occupied at T = K. One 
would hope that PP methods which treat all electrons 
self-consistently would be applicable for such finite-T cal- 
culations. A related issue is the validity of using PPs 
which remove some of the core. A rough estimate of the 
relevant scale comes from taking the Is ionization poten- 
tial for the Li atom to be approximately the magnitude of 
the LDA Kohn-Sham Is eigenvalue, about 51 eV. Then 
PP treatment of Li Is electrons as core might be expected 
to be applicable for temperatures much smaller than 51 
cV. The question is the validity of any estimate of this 
sort, in particular, how much smaller? We remark that 
Levashov et al. [36j found that for ambient density Al, 
the PAW pressure deviated from the all-electron value at 
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FIG. 6: Comparison of pressure vs. temperature for bcc-Li 
obtained with 3e— and le— pseudopotentials for the Perdew- 
Zunger LDA exchange-correlation functional as implemented 
in Quantum-Espresso (2-atom unit cell, 9x9x9 k-mesh, 128 
bands). Left panel: pLi = 0.5 g/cm 3 . Right panel: pLi = 1-0 
g/cm 3 . 



about T = 5-6 eV. This is less than 10% of the magnitude 
of the LSDA 2p atomic KS eigenvalue (about 70 eV) . 

First consider the comparative performance of the 
PPs. Figure [5] shows the hydrostatic pressure as a 
function of temperature calculated using le - and 3e _ 
norm-conserving pseudopotentials for the bcc-Li struc- 
ture (fixed nuclear positions, pu = 0.5 and 1.0 g/cm 3 ). 
Notice that this calculation uses a ground-state XC func- 
tional: there is no explicit temperature dependence in 
the PZ LDA XC functional. If the number of bands 
taken into account for a two-atom unit cell for pLi = 0.5 
g/cm -3 is 128, the occupation number of the highest 
energy bands is of the order of 10 -6 — 10 -7 . Observe 
that the results from the le - PP arc in almost perfect 
agreement with those from the 3e~ calculations for T 
up to 75,000 K, with small disagreement appearing at 
higher temperatures. For low to moderate compression, 
it appears that the range of applicability of standard le - 
norm-conserving pseudopotentials is at least up to T = 
100 kK or about 8-9 eV. This fits the rough argument 
based on the Li Is KS eigenvalue, with the criterion for 
"much smaller" being of order 20 % at most. 

Next comes the matter of significant fractional occupa- 
tion of ever-higher energy orbitals with increasing tem- 
perature. A related issue is energy level shifting and 
reordering with increased density, an effect known for 
T = OK Li dlSH. 

Again, we did calculations with ground-state XC on 
bcc Li, with material density from 0.6 to 4.0 g/cm 3 . For 
all temperatures and densities, we used a 7 x 7 x 7 
Monkhorst-Pack k-grid [50j . a two-atom unit cell, and 
included 128 bands with a plane wave energy cutoff of 
150 Ry. The calculations were done with Quantum- 
ESPRESSO and the le - PP just mentioned. The upper 
panel of Figure [7] shows a sample of the orbital eigen- 
values for a single k-point, T, as a function of material 
density for T = 100 kK. The T = 100 K plot is ab- 
solutely indistinguishable, due to relatively small differ- 
ences in the eigenvalues. The eigenvalues are labeled in 
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FIG. 7: Top: orbital energies for the T point, at 100 kK. 
Bottom left: single-spin occupation number fi of the levels 
plotted for T = 100 K. Bottom right: same but for T = 100 
kK. The legend for the level number, given in the upper right, 
is for all plots. 



order of increasing energy, e\ lowest, E\2b highest. The 
main point to be noticed is that as the density increases, 
the spread in the lower half (roughly) of the eigenval- 
ues increases. Those are the eigenvalues most pertinent 
to the calculation, in the sense that at a given tempera- 
ture, excited levels will be depopulated at higher densi- 
ties compared to the corresponding levels at lower densi- 
ties. (An exception would be a pressure-induced switch 
in level-ordering.) The lower panels of Fig. [7] show the 
occupation numbers for those same eigenvalues. At low 
temperature and low density, the results are as expected, 
an almost square-wave Fermi distribution with the low- 
est band fully occupied (since there are two electrons 
in the unit cell) and the higher bands unoccupied. At 
higher densities, some k-points, including the V point, as 
shown in the lower left panel, for densities above 3 g/ cm 3 , 
have no occupation while others have two occupied lev- 
els. This repopulation is a consequence of changes in 
the KS orbitals caused by changes in the external poten- 
tial, hence also in the effective KS potential. The lower 
right panel shows that at higher temperatures there is not 
only a temperature dependence of the occupation num- 
bers, but a significant density dependence because of the 
spreading of the orbital energy levels. 

Next we consider the number of bands required for 
a stipulated precision, given by a minimum occupation 
number threshold, as a function of temperature and den- 
sity. For the calculations just discussed, we calculated a 
zonc-avcraged band occupation. For a band of composite 
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FIG. 8: Zone-averaged band occupation numbers for all 
bands at T = 100 kK, for the various material densities listed. 



index i, we sum the occupations of the £j level multiplied 
by the k-point integration weight for all k-points. This 
zone-averaged occupation is plotted for T = lOOkK in 
Figure [8] One sees clearly that, for a given threshold in 
occupation number, for example 10 -6 , the number of re- 
quired bands decreases significantly with increasing den- 
sity. This decrease again is due to changes in the KS or- 
bitals. At least at T = OK, it long has been known 44 , 49[ 
that as Li is compressed, its band structure initially be- 
comes less like the homogeneous electron gas (HEG) than 
the bcc zero-pressure bands. However, eventually, the 
system passes over to a Thomas-Fcrmi-Dirac equation of 
state, signifying near-perfect but spread parabolic bands 
(see upper panel of Fig. [7]) and corresponding HEG oc- 
cupations. 



B. Exchange free energy 

Though it originated in the Greens function for- 
malism of many-fermion theory, the finite-tcmpcraturc 
Hartree-Fock approximation is the thermodynamical 
generalization of the variational optimization of a single- 
determinant trial wave function which is ubiquitous in 
quantum-chemistry and molecular physics as the HF ap- 
proximation (Hij . 

To summarize, the thermal generalization of the fa- 
miliar HF single-determinantal exchange energy may be 
expressed in terms of the one-electron reduced density 
matrix (1-RDM) 

•Fx[n] :=- J d Xl dx 2 {<? 12 f«(xi|x 2 ) 

xfW(x 2 |x' 1 )} xi=XliX - =X2 , (7) 

where x r, s is a composite space-spin variable, 512 = 
l/2|ri — r2 1 , and the 1-RDM is defined in terms of the 



relevant orbitals {tpi} and occupation numbers {fi} 

00 

f( 1 )(x 1 |x' 1 ):=E/i^'( x i)^K)> (8) 
i=i 

subject to 

= f(e j - M ) = [1 + exp(/3( £j - /i))]" 1 (9) 

and 

dxLpi(x.)ip*{x) = Sij 
N , 



(10) 



with f3 := 1/fc^T as usual. Here [i is the chemical poten- 
tial (determined by Eq. (fTt))0 and the Ej are the eigen- 
values of the associated one-particle ftHF equation. 

The analogue to ftHF in DFT is called finite- 
temperature exact exchange (ftEXX hereafter) DFT [52| . 
In its pure Kohn-Sham form, ftEXX defines the exchange 
free energy formally identically with ftHF, but evaluates 
the density n(r, T) from orbitals which follow from a true 
KS procedure, that is, from a one-body Hamiltonian with 
a local (multiplicative) exchange potential. That poten- 
tial follows from the system response function, dn/dvKS- 
A full ftDFT calculation (not exchange only) would have 
a correlation free energy functional and associated KS 
potential as well. 

Ground state DFT with so-called hybrid approximate 
exchange functionals has a similar structure for the total 
energy, in the sense that hybrids have contributions both 
from single-determinant exchange and from exchange- 
correlation functionals which are explicitly density de- 
pendent. Instead of a KS procedure, one can go from such 
a hybrid expression directly to coupled one-electron equa- 
tions by explicit variation with respect to the orbitals. In 
ground-state theory with a hybrid functional, this proce- 
dure sometimes is called generalized KS. The relevant 
point is that the same approach applies directly to ftHF. 
Simply switch off the explicit density functionals for ex- 
change and correlation and leave the exchange functional 
which comes from the trace over single-determinants. 
Since the capacity to do hybrid DFT as a generalized KS 
approach exists in both Vasp and Quantum-Espresso, 
one sees that such coding is immediately exploitable for 
doing ftHF. 

In parallel with ground state DFT, an LDA may be 
obtained from considering the finite-T HEG. Its exchange 
free energy is given in first-order perturbation theory by 



HEG 



V 



(27r) e 



47T 



dkdk'-— ^/(k)/(k') (11) 



where /(k) = [1 +exp(/3(fc 2 /2 - M HEG ))]-\ and V is the 
system volume. With the chemical potential expanded 
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to the same order as well, /i HEG 
exchange portion is 



/i™ G , the 



M X HEG KT) 



HEG 



(12) 



If expressed in closed form, this result may be used 
as the finite-T LDA, with exchange free energy per 
electron # DA (n(r, T), T) = (T^ EG /nV)\ 



,LDA/ 



n=n(r,T)> and 

Here we 



v-^(n(v,T),T) = 0™ G (n, T) | n=n(l ., T) . 
used the parametrization given by Perrot and Dharma- 
wardana [53[. The LDA exchange free energy is then 

7- x LDA Kr),T]= |/ x LDA (n(r,T),T)n(r,T)dr. (13) 

The one-particle density follows by obvious analogy with 

Eqs.u-rro] 



C. Finite-temperature Hartree-Fock and DFT 
X-only calculations 

To study the importance of using an explicitly T- 
dependcnt expression for the exchange free energy (rather 
than a calculation with a ground-state X functional) and 
to estimate the quality of the T-dependent exchange free- 
energy functional defined by Eq. (fT3|) . we compare ftHF 
calculations which use the exact exchange free energy Eq. 
([7|), Kohn-Sham calculations with T-independent LDA 
exchange for the exchange free energy, « E^ DA , and 
KS calculations done with the T-dependent exchange free 
energy functional J r x DA - In the following discussion, the 
ground-state functional calculations are labeled "LDAx" , 
while those which used the explicitly T-dependent LDA 
are labeled "LDAx(T)". All the calculations were done 
with Quantum-Espresso using the le~~ PZ LDA pscu- 
dopotential taken from the Quantum-Espresso web 
page. We treated bcc Li with fixed nuclear positions, here 
with densities between pLi = 0.6 and 1.8 g/cm 3 and tem- 
peratures between 100K and lOOkK. At these densities 
and temperatures, the Quantum-Espresso le~ pseu- 
dopotcntial is adequate; recall Sections IIII Bl and IIV Al as 
well as Figs. [5] and [6] 

Convergence of the ftHF and the LDAx calculations 
with respect to the k-mesh for the bcc-Li 2-atom unit 
cell requires attention. It is known (Hif that T=0 K LDA 
calculations on bcc Li exhibit misleading convergence be- 
havior at a relatively coarse k-mesh density. We tested 
for the smallest real space cell size used, corresponding to 
bulk density pu = 1.8 g/cm 3 . The ftHF total free energy 
calculations converge much more slowly than the ftDFT 
calculations. For the moderate 7x7x7 k-mesh, the DFT 
calculations are converged to an iteration-to-iteration dif- 
ference of 0.02 eV per atom, while the HF calculations 
converge to the same precision only upon reaching the 
much denser 17 x 17 x 17 mesh. Moreover, the HF calcu- 
lation exhibits a potentially misleading energy minimum 
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FIG. 9: Comparison of finite-temperature HF, ground 
state LDA ~K-only (LDAx) and T-dependent LDA yi-only 
(LDAx(T)) exchange free energy differences AJ r x (T) = 
J-" X (T)— J-x(lOOK) per atom as afunction of electronic temper- 
ature T. Left panel: pLi = 0.6 g/cm 3 ; right panel: pLi = 1.2 
g/cm 3 . 



at 15 x 15 x 15. The k-mesh convergence becomes faster 
with increasing T. For example, at T = 100 kK, both 
HF and LDAx calculations already are converged at the 
3x3x3 k-mesh. In all calculations, both HF and DFT, 
presented in this section, the 25 x 25 x 25 k-mesh was 
used. 

Figure [9] compares changes in the exchange free energy 
contribution with increasing T relative to 100K values, 
7" X (T) - J" X (100K). The T-independent LDA exchange 
free energy practically does not change over that range, 
i.e., j£ DA [n(r,T)] « j£ DA [7i(r, 100K)]. In contrast, the 
HF exchange free energy increases significantly (by about 
4-5 eV per atom) with increasing T. The T-dependent 
LDA exchange free energy reproduces the HF behavior 
at least qualitatively. 

The exchange free energy is, of course, a small portion 
of the total free energy. Figure [TBI shows total free energy 
differences AJ" tot (T) = J" tot (T) - J" tot (100K) as a func- 
tion of electronic temperature. The free energy is mono- 
tonically decreasing with increasing T, in agreement with 
non-negativity of the entropy evaluated from the thermo- 



dynamic relation S 



dT 



N,V 



The DFT X-only total 



free energies from T-indepcndcnt LDA X lie below the 
corresponding ftHF values for all T and both densities. 
At T « 20 kK, the interval is about 2 eV/atom, growing 
to about 4-5 eV/atom by 40 kK. The T-dependent LDA 
X gives total free energy behavior much closer to that of 
ftHF, with discrepancies not exceeding 1 — 2 eV/atom. 

The effect of explicit T-dependence in exchange upon 
the pressure may be estimated from the difference be- 
tween the ftHF or LDAx(T) and the LDAx values. Fig- 
ure [TTJ provides this comparison. As a function of T, the 
pressure from ftHF starts below the LDAx curve, then 
crosses and goes above it at about 55-75 kK, depend- 
ing upon the material density. The temperature of this 
crossing point increases slightly with material density. To 
isolate the effect of exact T-dcpcndcnt X upon the pres- 
sure, we consider the difference between ftHF and LDAx 
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FIG. 10: Comparison of finite-temperature HF, ground 
state LDA X-only (LDAx) and T-dependent LDA X- 
only (LDAx(T)) total free energy differences AJ r to t(T) = 
•7"tot (T) — J-tot (100K) per atom as a function of electronic tem- 
perature. Left panel: pLi = 0.6 g/cm 3 ; right panel: pu — 1.2 
g/cm 3 . 



values, offset by the near-zero-temperaturc difference 

APhf-ldAx(T) = Phf(T) - PldAx(T) (14) 

at T = 100K, i.e., 

AAPhf-ldAx(T) = APhf-ldAx(T) 

— APhF-LDAx (100K) . (15) 

One can see from the right-hand panels of Fig. QT] that 
the maximum magnitude of this difference at T= 100 
kK is about 10 % for all material densities considered. 
For low temperatures, the effect of exact T-dcpcndent 
X on pressure is stronger. For example, at 30 kK 
AAPnF-LDAx(30kK) m 5 GPa for material density 1.0 
g/cm 3 , the shift is ~ 30% of the HF pressure (about 
15 GPa) at that T. Again, the LDAx(T) and ftHF 
temperature-dependence resemble one another qualita- 
tively whereas the LDAx result does not. Note that the 
LDAx(T) crossing temperature with respect to the LDAx 
curve increases much more rapidly with increasing mate- 
rial density than for ftHF. For material density pLi =0.6 
g/cm 3 , both curves cross at T « 60 kK. At = 1.2 
g/cm 3 the LDAx(T) pressures crosses the LDAx curve 
at T k, 100 kK, higher than the temperature of the HF- 
LDAx crossing point T w 70 kK. With increasing ma- 
terial density the shift between LDAx(T) and ftHF in- 
creases especially for T > 50 kK. 



V. CONCLUSIONS 

Detailed computational examination of the applica- 
bility of standard PP and PAW methods to the WDM 
regime, with bulk Li as the test system, yields several 
insights. By unambiguous comparison with all-electron 
results from small Li clusters of bec-derived symmetry, 
we find that the PAW scheme requires a small augmen- 
tation sphere radius, that the compensation-charge term 
is not helpful, and that all electrons must be treated in 
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FIG. 11: Effect of temperature-dependent exchange on pres- 
sure. Left panels: pressure as a function of electronic tem- 
perature as predicted by HF, LDAx, and LDAx(T) calcu- 
lations for pLi = 0.6,0.8,1.0 and 1.2 g/cm 3 . Right pan- 
els: differences in pressure between calculations with In- 
dependent and T-independent X, P(HF) - P(LDAx) and 
P(LDAx(T)) -P(LDAx). 



the SCF calculation. We have constructed such PAW 
data sets for LDA and GGA functionals and used them 
to generate reference data. 

We have located the maximal material density of bulk 
bcc-Li usable for standard PPs in Vasp, Abinit, and 
Quantum-Espresso codes. And we have delineated the 
validity of using such PPs at high T by comparison of le~ 
and 3e~ PP results. The transferability of PPs and PAW 
data sets developed for near-equilibrium conditions to the 
WDM regime is conditional. At near-equilibrium densi- 
ties it appears to be acceptable, but not at high densities. 
Clearly, such transferability should not be assumed. 

With these issues settled, we have found that there is 
non-trivial effect of explicit T-dcpcndence in the X func- 
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tional in the specific sense of comparison with ftHF. In 
particular, the LDA T-depcndcnt exchange contribution 
to the total free energy is much closer to the exact HF 
exchange value than is the contribution from exchange 
approximated by the LDA ground-state X functional. Al- 
though the exchange free energy is a small portion of the 
total free energy, this difference carries over into clearly 
significant differences in the equation of state. Thus, the 
effect of explicit T-dependence in X is relevant for an 
accurate characterization of the Li equation of state in 
the WDM regime. We suspect that this may be gen- 
erally true of WDM systems. If so, T-dependent LDA 
exchange may serve as a starting point for development 



of more refined GGA-type exchange free energy function- 
al, analogous with the role of LDA in the ground state. 
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